/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

/*******************************************************************************************************************************
 * sysLoad.c -- This module collects various functions employed to load an existing configuration.
 *******************************************************************************************************************************/

#include "dataStructs.h"

extern VecR *siteCoords, *siteOrientVec, *siteVelocity, *atomCoordinates, *siteAngMom; 
extern RMat *siteRotMat;
extern int nPointMasses, nSymmRigBodies, nNonSymmRigBodies;
extern const int zConstraint, nSolutes;
extern Atom *atom;  /* array of Atom, represents solute atoms */
extern Site *site;  /* array of Site, represents the state of the system */

/* external function declarations */
void EulerToRotMat( RMat*, real, real, real );
void GetOrientPrincAxisZ_FromRotMat( VecR*, RMat );
real RandR();

void CheckStructureLoaded( const int nSites, const int nWaters ) // check the checkpoint file
{ 
  int n, nPointMassSites, nSymmRigBodySites, nWats, nChol = 0, nAmine = 0, nPhos = 0, nGly = 0, nEst = 0, nTail = 0, nWat = 0, nSolute = 0, totNSites;
  DO_SITE {
    switch ( site[ n ].lipidUnit ) {
    case 0: 
      nWat++; 
      break;
    case 1: 
      if ( site[ n ].type == CHOLINE_TYPE )
	nChol++; 
      else if ( site[ n ].type == AMINE_TYPE )
	nAmine++; 
      break;
    case 2: 
      nPhos++; 
      break;
    case 3: 
      nGly++; 
      break;
    case 7: 
      nGly++; 
      break;
    case 4: 
    case 10:
      nEst++;
      break;
    case 5: 
    case 6: 
    case 8: 
    case 9: 
    case 11:
    case 12:
    case 13:
    case 14:
    case 15:
      nTail++; 
      break;
    case SOLUTE: 
      nSolute++; 
      printf( "Loading solute # %d...\n", nSolute ); break;
    default: printf( "*** ERROR in CheckStructureLoaded(): got to default case.\n" ); exit(0);
    }
  }
  printf( "System (from checkpoint): %2d water, %2d choline, %2d amine, %2d phosphate, %2d glycerol, %2d ester and %2d tail sites ", 
	  nWat, nChol, nAmine, nPhos, nGly, nEst, nTail );
  if ( zConstraint && nSolutes ) printf( "and %d solute", nSolute );
  printf( "\n" );
  totNSites = nWat + nChol + nAmine + nPhos + nGly + nEst + nTail + nSolute;
  nPointMassSites = nChol + nAmine + nPhos + nTail; //printf( "#pointMasses: %d\t", nPointMassSites );
  nSymmRigBodySites = nGly + nEst + nWat; //printf( "#symmetricBodies: %d\t", nSymmRigBodySites );
  nWats = nWat; //printf( "#nonSymmBodies: %d\t\t", nWats + nSolute );
  if ( totNSites != nSites || totNSites != nPointMassSites + nSymmRigBodySites + nSolute || nWats != nWaters ) { 
    printf("\nWarning: #sites loaded (%d) != nSites (%d) (check your inputs!)\n", totNSites, nSites ); 
    if ( !zConstraint ) exit( 0 ); 
  };
  printf( "Total number of sites: %3d\n", totNSites );
}

void LoadSoluteParameters( const int insertSolute, const int nAtoms )
{
  FILE *fPtr;
  char line[ 81 ], *token;
  double d;
  int a;
  
  if ( ( fPtr = fopen( "solute.in", "r" ) ) == 0 ) {
    printf( "*** error: Can't open file solute.in!\nExiting...\n" );
    exit( 0 );
  }
  DO_ATOM {
    fgets( line, 81, fPtr );
    token = strtok( line, " \t\n" ); // find first token 
    strcpy( atom[ a ].type, token );
    token = strtok( NULL, " \t\n"  ); // 
    d = atof( token );
    atomCoordinates[ a ].x = d;
    token = strtok( NULL, " \t\n"  ); // 
    d = atof( token );
    atomCoordinates[ a ].y = d;
    token = strtok( NULL, " \t\n"  ); // 
    d = atof( token );
    atomCoordinates[ a ].z = d;
    token = strtok( NULL, " \t\n"  ); // 
    d = atof( token );
    atom[ a ].charge = d;
    token = strtok( NULL, " \t\n"  ); // 
    d = atof( token );
    atom[ a ].sig = d;
    token = strtok( NULL, " \t\n"  ); // 
    d = atof( token );
    atom[ a ].eps = d;
  }
  fclose( fPtr );

  if ( insertSolute ) {
    printf("\n Inserting & growing solute...\n");
    DO_ATOM {
      VScale( atomCoordinates[ a ], 0 );
      atom[ a ].charge *= 0;
      atom[ a ].sig *= 0;
      atom[ a ].eps *= 0;
    }
  }
  printf("\nSolute loaded:\n");
  printf( "Atom\t\tCoords/A\t Charge/e  SigLJ/A  EpsLJ/kcal/mol\n" );
  DO_ATOM 
    printf( "%5s\t%7.3f%7.3f%7.3f\t%9.5f%9.5f%9.5f\n", atom[a].type, atomCoordinates[a].x, atomCoordinates[a].y, atomCoordinates[a].z, atom[a].charge, atom[a].sig, atom[a].eps );      
}

void LoadVelFile( const int nSites ) {
  FILE *fp;
  int lipId, i = 0;
  //nChol = 0, nPhos = 0, nGly = 0, nTail = 0, nWat = 0, n, totNSites;
  char line[ 80 ], *tokenPtr;
  printf("*** loading velocity file...\n");
  if ( ( fp = fopen( "velocities.out", "r" ) ) == 0 ) printf( "*** error: Can't open velocity file!\n" );
  else{
    while ( i < nSites ) {
      fgets( line, 80, fp ); /* read one line and store it in "line" */
      tokenPtr = strtok( line, " \t\n" ); /* break the string "line" into the substrings delimited by any of " \t\n"  */
      lipId = atoi ( tokenPtr );
      
      switch ( site[ i ].mechType ) {
      case PM: 
	tokenPtr = strtok( NULL, " \t\n" );
	siteVelocity[ i ].x = atof( tokenPtr ); /* load COM velocities */
	tokenPtr = strtok( NULL, " \t\n" );
	siteVelocity[ i ].y = atof( tokenPtr ); /* load COM velocities */
	tokenPtr = strtok( NULL, " \t\n" );
	siteVelocity[ i ].z = atof( tokenPtr ); /* load COM velocities */
	break;      
      case SRB: /* fall through */
      case NSRB:
	tokenPtr = strtok( NULL, " \t\n" );
	siteVelocity[ i ].x = atof( tokenPtr ); /* load COM velocities */
	tokenPtr = strtok( NULL, " \t\n" );
	siteVelocity[ i ].y = atof( tokenPtr ); /* load COM velocities */
	tokenPtr = strtok( NULL, " \t\n" );
	siteVelocity[ i ].z = atof( tokenPtr ); /* load COM velocities */
	tokenPtr = strtok( NULL, " \t\n" );
	siteAngMom[ i ].x = atof( tokenPtr ); /* load angular momenta */
	tokenPtr = strtok( NULL, " \t\n" );
	siteAngMom[ i ].y = atof( tokenPtr ); /* load angular momenta */
	tokenPtr = strtok( NULL, " \t\n" );
	siteAngMom[ i ].z = atof( tokenPtr ); /* load angular momenta */
	break;
      default: printf( "*** ERROR in LoadVelFile(): got to default case...\n" );
      }
      if ( tokenPtr == NULL ) break;
      ++i;
    } /* end while */
    fclose( fp );
  }
}

void LoadStructureFile( const int nSites, const int nWaters )
{
  FILE *fp;
  int i = 0, nChol = 0, nPhos = 0, nGly = 0, nTail = 0, nWat = 0, n, nWats, totNSites;
  char line[ 80 ], *tokenPtr;
  real eAng[ 3 ], ez, ey;
  /*   fPtr = fopen("structure.log","w"); */
  printf("*** loading structure file...\n");
  if ( ( fp = fopen( "structure.in", "r" ) ) == 0 ) printf( "*** error: Can't open structure file!\n" );
  else{
    while ( i < nSites ) {
      fgets( line, 80, fp ); /* read one line and store it in "line" */
      /* fprintf( fPtr, "%s", line ); */
      /*printf("%s", line);*/
      /*       if ( feof( fp ) ) break;  */
      tokenPtr = strtok( line, " \t\n" ); /* break the string "line" into the substrings delimited by any of " \t\n"  */
      site[ i ].lipidUnit = atoi ( tokenPtr );
      switch ( site[ i ].lipidUnit ) {
      case 0: site[ i ].type = WATER_TYPE; site[ i ].mechType = NSRB; nWat++; break;
      case 1: site[ i ].type = CHOLINE_TYPE; site[ i ].mechType = PM; nChol++; break;
      case 2: site[ i ].type = PHOSPHATE_TYPE; site[ i ].mechType = PM; nPhos++; break;
      case 3: /* fall through */
      case 4: /* fall through */
      case 5: /* fall through */
      case 6: /* fall through */
      case 8: /* fall through */
      case 9: /* fall through */
      case 10:
      case 11:
      case 12: site[ i ].type = TAIL_TYPE; site[ i ].mechType = SRB; nTail++; break;
      default: printf( "*** ERROR in LoadStructureFile(): got to default case...\n" ); exit(0);
      }
      tokenPtr = strtok( NULL, " \t\n" );
      siteCoords[ i ].x = atof( tokenPtr ); /* load COM coordinates */
      tokenPtr = strtok( NULL, " \t\n" );
      siteCoords[ i ].y = atof( tokenPtr ); /* load COM coordinates */
      tokenPtr = strtok( NULL, " \t\n" );
      siteCoords[ i ].z = atof( tokenPtr ); /* load COM coordinates */
      /*     if ( site[ i ].mechType == 1 ) {/\* if SRB lets also load the orientation vector *\/ */
      if ( site[ i ].mechType == SRB ) {/* if SRB lets also load the orientation vector */
	tokenPtr = strtok( NULL, " \t\n" );
	siteOrientVec[ i ].x = atof( tokenPtr ); /* load SRB orientation */
	tokenPtr = strtok( NULL, " \t\n" );
	siteOrientVec[ i ].y = atof( tokenPtr ); /* load SRB orientation */
	tokenPtr = strtok( NULL, " \t\n" );
	siteOrientVec[ i ].z = atof( tokenPtr ); /* load SRB orientation */

	//	printf("*** e loaded   = % f, % f, % f\n", siteOrientVec[ i ].x, siteOrientVec[ i ].y, siteOrientVec[ i ].z );

	//	GetRotMatFromOrientVector( &siteRotMat[ i ], &siteOrientVec[ i ] ); 

	//	GetOrientPrincAxisZ_FromRotMat( &siteOrientVec[ i ], siteRotMat[ i ] );
	
	// trick to get the correct sign in e_z
	ez = siteOrientVec[ i ].z;
	ey = siteOrientVec[ i ].y;

	eAng[ 1 ] = - atan2( siteOrientVec[ i ].x, siteOrientVec[ i ].y );
	eAng[ 2 ] = - acos( siteOrientVec[ i ].z );
	eAng[ 0 ] = 2. * PI * RandR();

	/* 	EulerToQuat( &qe, eAng ); */
	/* 	BuildRotMatrixFromQuat( &siteRotMat[ i ], &qe, 1 ); */

	EulerToRotMat( &siteRotMat[ i ], eAng[1], eAng[2], eAng[0] );
	GetOrientPrincAxisZ_FromRotMat( &siteOrientVec[ i ], siteRotMat[ i ] ); /* convert rotation matrix to orientation vector */

	if ( ez * siteOrientVec[ i ].z < 0 || ey * siteOrientVec[ i ].y < 0 || ez * siteOrientVec[ i ].z < 0 ) { // if the original sign and the converted dont match
	  printf("WARNING: orientation in structure.in not converted properly to the actual orientation!\n");
	  /* 	  siteOrientVec[ i ].z = - siteOrientVec[ i ].z; */
	}
	/* 	if ( ey * siteOrientVec[ i ].y < 0 ) { // if the original sign and the converted dont match */
	/* 	  siteOrientVec[ i ].y = - siteOrientVec[ i ].y; */
	/* 	} */

	//	printf("*** e converted = % f, % f, % f\n\n", siteOrientVec[ i ].x, siteOrientVec[ i ].y, siteOrientVec[ i ].z );

      }
      if ( tokenPtr == NULL ) break;
      ++i;
    } /* end while */

    fclose( fp );
  }

  printf( "Structure loaded: %2d waters, %2d cholines, %2d phosphates, %2d glycerols and %2d tail sites\n", nWat, nChol, nPhos, nGly, nTail );
  totNSites = nWat + nChol + nPhos + nGly + nTail;
  nPointMasses = nChol + nPhos; printf( "#pointMasses: %d\t", nPointMasses );
  nSymmRigBodies = nGly + nTail; printf( "#symmetricBodies: %d\t", nSymmRigBodies );
  nWats = nWat; printf( "#nonSymmBodies: %d\t\t", nWats );
  if ( totNSites != nSites || totNSites != nPointMasses + nSymmRigBodies + nWats || nWats != nWaters ) { 
    printf("\nWarning: #sites loaded != nSites (check your inputs!)\n"); 
    if ( !zConstraint || !nSolutes ) exit( 0 ); 
  };
  DO_SITE VScale( siteCoords[ n ], Angstrom_IN_BRAHMS_UNITS ); /* scaling coordinates to reduced units */
}


